# Load packages
  library(tidyverse)
  library(broom)
  library(lmtest)
  library(sandwich)
  library(stargazer)

# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]

  #### MEDW 2016 Spanish National Election Study ----
  # Load data
  load("ObsStudy_00_data_Spain_2016.RData")
  dt <- table

  # Party rating
  colnames(dt)[seq(47, 55, 1)] <- c("rat_PP", "rat_PSOE", "rat_Podemos", "rat_Ciudadanos", "rat_IU", "rat_CDC", "rat_ERC", "rat_EAJ-PNV", "rat_EHBildu")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 47:55][dt[, 47:55] == 99] <- NA
  
  dt$rat_PP <- as.numeric(dt$rat_PP)
  dt$rat_PSOE <- as.numeric(dt$rat_PSOE)
  dt$rat_Podemos <- as.numeric(dt$rat_Podemos)
  dt$rat_Ciudadanos <- as.numeric(dt$rat_Ciudadanos)
  
  # Coalition evaluations/preferred coalition
  colnames(dt)[seq(70, 74, 1)] <- c("rat_coal_PP_PSOE", "rat_coal_PP_Ciudadanos", "rat_coal_PSOE_Podemos", "rat_coal_PSOE_Ciudadanos", "rat_coal_Podemos_Ciudadanos")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 70:74][dt[, 70:74] == 99] <- NA
  
  dt$rat_coal_PP_PSOE <- as.numeric(dt$rat_coal_PP_PSOE)
  dt$rat_coal_PP_Ciudadanos <- as.numeric(dt$rat_coal_PP_Ciudadanos)
  dt$rat_coal_PSOE_Podemos <- as.numeric(dt$rat_coal_PSOE_Podemos)
  dt$rat_coal_PSOE_Ciudadanos <- as.numeric(dt$rat_coal_PSOE_Ciudadanos)
  dt$rat_coal_Podemos_Ciudadanos <- as.numeric(dt$rat_coal_Podemos_Ciudadanos)
    
  
  
  # Coalition likelihood
  colnames(dt)[seq(76, 80, 1)] <- c("coallik_PP_PSOE", "coallik_PP_Ciudadanos", "coallik_PSOE_Podemos", "coallik_PSOE_Ciudadanos", "coallik_Podemos_Ciudadanos")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 76:80][dt[, 76:80] == 99] <- NA
  
  # Dependent Variable: nominal vote choice
  # Set missings (Other Party (88), Blank Vote (89), Null Vote (90), Don´t know (99)) to NA 
  dt[, 25][dt[, 25] >= 88] <- NA
  
  # Vote variable
  dt$vote <- dt$Q8A
  
  dt$vote[dt$vote==1] <- "PP"
  dt$vote[dt$vote==2] <- "PSOE"
  dt$vote[dt$vote==3] <- "Podemos"
  dt$vote[dt$vote==4] <- "Ciudadanos"
  dt$vote[dt$vote==5] <- "ECP"
  dt$vote[dt$vote==6] <- "EUPV"
  dt$vote[dt$vote==7] <- "ERC"
  dt$vote[dt$vote==8] <- "CDC"
  dt$vote[dt$vote==9] <- "EnMarea"
  dt$vote[dt$vote==10] <- "EAJ-PNV"
  dt$vote[dt$vote==11] <- "EHBildu"
  dt$vote[dt$vote==12] <- "CC"
  dt$vote[dt$vote==13] <- "UPN"
  dt$vote[dt$vote==14] <- "PAR"
  dt$vote[dt$vote==15] <- "ForoAsturias"
  dt$vote[dt$vote==16] <- "CHA"
  dt$vote[dt$vote==17] <- "UPyD"
  dt$vote[dt$vote==18] <- "BNG"
  dt$vote[dt$vote==19] <- "GeroaBai"
  dt$vote[dt$vote==20] <- "PACMA"
  
  # Generate ptv columns based on vote variable
  dt$ptv_PP <- ifelse(dt$vote == "PP", 1, 0)
  dt$ptv_PSOE <- ifelse(dt$vote == "PSOE", 1, 0)
  dt$ptv_Podemos <- ifelse(dt$vote == "Podemos", 1, 0)
  dt$ptv_Ciudadanos <- ifelse(dt$vote == "Ciudadanos", 1, 0)
  
  # Coalition evaluation
  class(dt$rat_coal_PP_PSOE)
  dt$rat_coal_PP_PSOE <- as.numeric(dt$rat_coal_PP_PSOE)
  
  class(dt$rat_coal_PP_Ciudadanos)
  dt$rat_coal_PP_Ciudadanos <- as.numeric(dt$rat_coal_PP_Ciudadanos)
  
  class(dt$rat_coal_PSOE_Podemos)
  dt$rat_coal_PSOE_Podemos <- as.numeric(dt$rat_coal_PSOE_Podemos)
  
  class(dt$rat_coal_PSOE_Ciudadanos)
  dt$rat_coal_PSOE_Ciudadanos <- as.numeric(dt$rat_coal_PSOE_Ciudadanos)
  
  class(dt$rat_coal_Podemos_Ciudadanos)
  dt$rat_coal_Podemos_Ciudadanos <- as.numeric(dt$rat_coal_Podemos_Ciudadanos)
  
  Z <- dt %>% select(contains("rat_coal")) %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_PP <- Z %>% select(rat_coal_PP_PSOE, rat_coal_PP_Ciudadanos) %>% as.matrix()
  rat_coal_PSOE <- Z %>% select(rat_coal_PP_PSOE, rat_coal_PSOE_Podemos, rat_coal_PSOE_Ciudadanos) %>% as.matrix()
  rat_coal_Podemos <- Z %>% select(rat_coal_PSOE_Podemos, rat_coal_Podemos_Ciudadanos) %>% as.matrix()
  rat_coal_Ciudadanos <- Z %>% select(rat_coal_PP_Ciudadanos, rat_coal_PSOE_Ciudadanos, rat_coal_Podemos_Ciudadanos) %>% as.matrix()
  
  
  # Coalition likelihood
  gamma <- dt %>% select(contains("coallik_")) 
  
  coal_lik_PP <- gamma %>% select(coallik_PP_PSOE, coallik_PP_Ciudadanos) %>%
    mutate(sum_exp = exp(coallik_PP_PSOE) + exp(coallik_PP_Ciudadanos),
           coallik_PP_PSOE = exp(coallik_PP_PSOE)/sum_exp,
           coallik_PP_Ciudadanos = exp(coallik_PP_Ciudadanos)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_PSOE <- gamma %>% select(coallik_PP_PSOE, coallik_PSOE_Podemos, coallik_PSOE_Ciudadanos) %>%
    mutate(sum_exp = exp(coallik_PP_PSOE) + exp(coallik_PSOE_Podemos) + exp(coallik_PSOE_Ciudadanos),
           coallik_PP_PSOE = exp(coallik_PP_PSOE)/sum_exp,
           coallik_PSOE_Podemos = exp(coallik_PSOE_Podemos)/sum_exp,
           coallik_PSOE_Ciudadanos = exp(coallik_PSOE_Ciudadanos)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_Podemos <- gamma %>% select(coallik_PSOE_Podemos, coallik_Podemos_Ciudadanos) %>%
    mutate(sum_exp = exp(coallik_PSOE_Podemos) + exp(coallik_Podemos_Ciudadanos),
           coallik_PSOE_Podemos = exp(coallik_PSOE_Podemos)/sum_exp,
           coallik_Podemos_Ciudadanos = exp(coallik_Podemos_Ciudadanos)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_Ciudadanos <- gamma %>% select(coallik_PP_Ciudadanos, coallik_PSOE_Ciudadanos, coallik_Podemos_Ciudadanos) %>%
    mutate(sum_exp = exp(coallik_PP_Ciudadanos) + exp(coallik_PSOE_Ciudadanos) + exp(coallik_Podemos_Ciudadanos),
           coallik_PP_Ciudadanos = exp(coallik_PP_Ciudadanos)/sum_exp,
           coallik_PSOE_Ciudadanos = exp(coallik_PSOE_Ciudadanos)/sum_exp,
           coallik_Podemos_Ciudadanos = exp(coallik_Podemos_Ciudadanos)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  # Mean-variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  
  M_PP <- sapply(1:N, function(i) EV(coal_lik_PP[i,], rat_coal_PP[i,]))
  V_PP <- sapply(1:N, function(i) Vari(coal_lik_PP[i,], rat_coal_PP[i,]))
  
  M_PSOE <- sapply(1:N, function(i) EV(coal_lik_PSOE[i,], rat_coal_PSOE[i,]))
  V_PSOE <- sapply(1:N, function(i) Vari(coal_lik_PSOE[i,], rat_coal_PSOE[i,]))
  
  M_Podemos <- sapply(1:N, function(i) EV(coal_lik_Podemos[i,], rat_coal_Podemos[i,]))
  V_Podemos <- sapply(1:N, function(i) Vari(coal_lik_Podemos[i,], rat_coal_Podemos[i,]))
  
  M_Ciudadanos <- sapply(1:N, function(i) EV(coal_lik_Ciudadanos[i,], rat_coal_Ciudadanos[i,]))
  V_Ciudadanos <- sapply(1:N, function(i) Vari(coal_lik_Ciudadanos[i,], rat_coal_Ciudadanos[i,])) 
  
  
  # Gender
  dt$male <- ifelse(dt$GEND == 1, 1, 0)

  # Age
  dt$AGE <- as.numeric(dt$AGE)
  
  # Education
  dt$SD4 <- as.numeric(dt$SD4)
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$pid <- dt$Q47A
  dt$pid[dt$Q47==2] <- 0
  
  dt$pid_PP <- ifelse(dt$pid==1, 1, 0)
  dt$pid_PSOE <- ifelse(dt$pid==2, 1, 0)
  dt$pid_Podemos <- ifelse(dt$pid==3, 1, 0)
  dt$pid_Ciudadanos <- ifelse(dt$pid==4, 1, 0)
  
  d <- data.frame(
    "ptv_PP" = dt$ptv_PP,
    "ptv_PSOE" = dt$ptv_PSOE,
    "ptv_Podemos" = dt$ptv_Podemos,
    "ptv_Ciudadanos" = dt$ptv_Ciudadanos,
    "rat.PP" = dt$rat_PP,
    "rat.PSOE" = dt$rat_PSOE,
    "rat.Podemos" = dt$rat_Podemos,
    "rat.Ciudadanos" = dt$rat_Ciudadanos,
    "pid_PP" = dt$pid_PP,
    "pid_PSOE" = dt$pid_PSOE,
    "pid_Podemos" = dt$pid_Podemos,
    "pid_Ciudadanos" = dt$pid_Ciudadanos,
    "lotterymean.PP" = M_PP,
    "lotteryvariance.PP" = V_PP,
    "lotterymean.PSOE" = M_PSOE,
    "lotteryvariance.PSOE" = V_PSOE,
    "lotterymean.Podemos" = M_Podemos,
    "lotteryvariance.Podemos" = V_Podemos,
    "lotterymean.Ciudadanos" = M_Ciudadanos,
    "lotteryvariance.Ciudadanos" = V_Ciudadanos,
    "sex" = dt$male,
    "age" = dt$AGE,
    "age_sqrt" = (dt$AGE)^2,
    "edu" = dt$SD4,
    "rat_coal_PP_PSOE" = dt$rat_coal_PP_PSOE,
    "rat_coal_PP_Ciudadanos" = dt$rat_coal_PP_Ciudadanos,
    "rat_coal_PSOE_Podemos" = dt$rat_coal_PSOE_Podemos,
    "rat_coal_PSOE_Ciudadanos" = dt$rat_coal_PSOE_Ciudadanos,
    "rat_coal_Podemos_Ciudadanos" = dt$rat_coal_Podemos_Ciudadanos
  )

  d$lotteryvariance.PP <- scales::rescale(d$lotteryvariance.PP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.PSOE <- scales::rescale(d$lotteryvariance.PSOE, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.Podemos <- scales::rescale(d$lotteryvariance.Podemos, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.Ciudadanos <- scales::rescale(d$lotteryvariance.Ciudadanos, to = c(0, 1), from = c(0,0.25))
  
# Save model results
  summary(m1 <- lm(as.formula(paste("ptv_PP ~ lotteryvariance.PP + lotterymean.PP + rat.PP + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_PP" else "")), d))
  summary(m2 <- lm(as.formula(paste("ptv_PSOE ~ lotteryvariance.PSOE + lotterymean.PSOE + rat.PSOE + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_PSOE" else "")), d))
  summary(m3 <- lm(as.formula(paste("ptv_Podemos ~ lotteryvariance.Podemos + lotterymean.Podemos + rat.Podemos + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_Podemos" else "")), d))
  summary(m4 <- lm(as.formula(paste("ptv_Ciudadanos ~ lotteryvariance.Ciudadanos + lotterymean.Ciudadanos + rat.Ciudadanos + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_Ciudadanos" else "")), d))
  
  ESP_2016_SNES_PP_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.PP") %>% select("estimate")
  ESP_2016_SNES_PP_Variance_SE <- se_robust(m1)["lotteryvariance.PP"] #tidy(m1) %>% filter(term=="lotteryvariance.PP") %>% select("std.error")
  ESP_2016_SNES_PP_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.PP") %>% select("estimate")
  ESP_2016_SNES_PP_Mean_SE <- se_robust(m1)["lotterymean.PP"] #tidy(m1) %>% filter(term=="lotterymean.PP") %>% select("std.error")
  
  ESP_2016_SNES_PSOE_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.PSOE") %>% select("estimate")
  ESP_2016_SNES_PSOE_Variance_SE <- se_robust(m2)["lotteryvariance.PSOE"] #tidy(m2) %>% filter(term=="lotteryvariance.PSOE") %>% select("std.error")
  ESP_2016_SNES_PSOE_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.PSOE") %>% select("estimate")
  ESP_2016_SNES_PSOE_Mean_SE <- se_robust(m2)["lotterymean.PSOE"] #tidy(m2) %>% filter(term=="lotterymean.PSOE") %>% select("std.error")
  
  ESP_2016_SNES_Podemos_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.Podemos") %>% select("estimate")
  ESP_2016_SNES_Podemos_Variance_SE <- se_robust(m3)["lotteryvariance.Podemos"] #tidy(m3) %>% filter(term=="lotteryvariance.Podemos") %>% select("std.error")
  ESP_2016_SNES_Podemos_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.Podemos") %>% select("estimate")
  ESP_2016_SNES_Podemos_Mean_SE <- se_robust(m3)["lotterymean.Podemos"] #tidy(m3) %>% filter(term=="lotterymean.Podemos") %>% select("std.error")
  
  ESP_2016_SNES_Ciudadanos_Variance_Estimate <- tidy(m4) %>% filter(term=="lotteryvariance.Ciudadanos") %>% select("estimate")
  ESP_2016_SNES_Ciudadanos_Variance_SE <- se_robust(m4)["lotteryvariance.Ciudadanos"] #tidy(m4) %>% filter(term=="lotteryvariance.Ciudadanos") %>% select("std.error")
  ESP_2016_SNES_Ciudadanos_Mean_Estimate <- tidy(m4) %>% filter(term=="lotterymean.Ciudadanos") %>% select("estimate")
  ESP_2016_SNES_Ciudadanos_Mean_SE <- se_robust(m4)["lotterymean.Ciudadanos"] #tidy(m4) %>% filter(term=="lotterymean.Ciudadanos") %>% select("std.error")
  
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.PP"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.PP"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.PSOE"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.PSOE"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.Podemos"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.Podemos"] <- "Government Lottery Mean"
  names(m4$coefficients)[names(m4$coefficients) == "lotteryvariance.Ciudadanos"] <- "Government Lottery Variance"
  names(m4$coefficients)[names(m4$coefficients) == "lotterymean.Ciudadanos"] <- "Government Lottery Mean"
  
  
  if(!robustnesscheck_pid){
    stargazer(m1, m2, m3, m4, title="Linear regressions of vote choice on perceived government lottery variance and mean (Spanish National Election Study 2016)", 
              align=TRUE, 
              dep.var.labels=c("PP", "PSOE", "Podemos", "Ciudadanos"),
              omit.stat=c("LL","ser","f"), type = "text",
              se = lapply(list(m1, m2, m3, m4), se_robust),
              p = lapply(list(m1, m2, m3, m4), p_robust),
              out = "TableSM7.tex"
    )
  }
 
  
